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1  INTRODUCTION 


1  Introduction 

We  may  say  that  one  of  the  main  objectives  of  the  kinetic  theory  is  to  describe  the  macro¬ 
scopic  properties  of  gases  -  like  pressure,  temperature,  thermal  conductivity,  viscosity, 
diffusion,  etc.  -  from  microscopic  quantities  that  are  associated  with  the  molecules  which 
compose  the  gases  -  like  mass,  velocity,  kinetic  energy,  internal  degrees  of  freedom  and 
interaction  forces  between  the  molecules. 

The  foundations  of  the  modern  kinetic  theory  of  gases  were  established  in  1867  by 
James  Clerk  Maxwell  (1831-1879)  who  proposed  a  general  transport  equation  for  arbi¬ 
trary  macroscopic  quantities  associated  with  mean  values  of  microscopic  quantities.  This 
equation  of  transport  relates  the  time  evolution  of  a  macroscopic  quantity  with  the  motion 
of  the  molecules,  collision  between  the  molecules  and  action  of  external  forces.  Although 
the  theory  was  valid  for  any  molecular  interaction  potential,  Maxwell  could  only  deter¬ 
mine  the  transport  coefficients  of  shear  viscosity,  thermal  conductivity  and  diffusion  by 
assuming  that  the  interaction  potential  was  derived  from  a  repulsive  force  which  was  in¬ 
versely  proportional  to  the  fifth  power  of  the  relative  distance  between  the  molecules. 
Nowadays  this  type  of  potential  is  known  as  Maxwellian  potential.  The  kinetic  theory  of 
gases  gained  a  new  impulse  in  1872  with  the  work  by  Ludwig  Eduard  Boltzmann  (1844- 
1906),  who  proposed  an  integro-differential  equation  -  the  Boltzmann  equation  -  which 
represents  the  evolution  of  the  velocity  distribution  function  in  the  phase  space  spanned 
by  the  coordinates  and  velocities  of  the  molecules.  In  the  Boltzmann  equation  the  tempo¬ 
ral  change  of  the  distribution  function  has  two  terms,  one  of  them  is  a  drift  term  due  to 
the  motion  of  the  molecules  while  the  other  one  is  a  collision  term  related  to  encounters 
of  the  molecules.  Based  on  this  equation,  Boltzmann  proposed  the  so-called  T~L  function 
which  decreases  with  time  or  remains  constant.  The  identification  of  this  function  as  the 
negative  of  the  gas  entropy  gave  a  molecular  interpretation  of  the  increase  of  the  entropy 
for  irreversible  processes.  Furthermore,  Boltzmann  in  the  same  work  presented  a  rigorous 
deduction  of  the  Maxwellian  distribution  function. 

From  the  Boltzmann  equation  one  could  determine  the  velocity  distribution  function 
hence  the  transport  coefficients  of  rarefied  gases,  however  this  task  was  not  so  easy.  Indeed, 
it  took  almost  forty  years  after  the  proposition  of  Boltzmann  equation,  for  David  Hilbert 
(1862-1943)  to  show  how  one  could  get  an  approximate  solution  of  the  integro-differential 
equation  from  a  power  series  expansion  of  a  parameter  which  is  proportional  to  the  mean 
free  path.  Further  advances  were  due  to  Sydney  Chapman  (1888-1970)  and  David  Enskog 
(1884-1947)  who  -  in  the  years  1916  and  1917  -  calculated  independently  and  by  different 
methods  the  transport  coefficients  for  gases  whose  molecules  interact  according  to  any  kind 
of  spherically  symmetric  potential  function.  Another  method  derived  from  the  Boltzmann 
equation  was  proposed  in  1949  by  Harold  Grad  (1923-1986)  who  expanded  the  distribution 
function  in  terms  of  tensorial  Hermite  polynomials  and  introduced  balance  equations 
corresponding  to  higher  order  moments  of  the  distribution  function. 

The  aim  of  these  notes  is  to  discuss  the  methods  of  Chapman-Enskog  and  Grad  with 
an  applications  to  granular  gases. 

(a)  Cartesian  notation  for  tensors  is  used  with  Latin  indexes  i,j,  k, . . .  -  which  may  range 
from  1  to  3  -  denoting  the  three-dimensional  system  of  spatial  coordinates  x,  y,  z; 

(b)  Einstein’s  summation  convention  over  repeated  indexes  is  used,  for  example,  T^Vj  = 
zCj=l  Tijvj'i 
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(c)  Parentheses  around  the  indexes  denote  the  symmetric  part  of  a  tensor,  brackets  its 
antisymmetric  part  while  angular  parentheses  refer  to  its  traceless  symmetric  part,  for 
example, 

+  Tji ), 

Above,  dij  is  Kronecker’s  symbol  and  the  traceless  tensor  is  called  tensor  deviator. 


i]  2 


T, 


T(ij)  T( 


(ij) 


It  a 

2  rr®ij  • 


■^(o)  2 


2  Boltzmann  and  Balance  Equations 


We  consider  a  monatomic  rarefied  gas  with  N  molecules  enclosed  in  a  recipient  of  volume 
V,  where  a  molecule  is  specified  as  a  point  in  a  six-dimensional  phase  space  spanned  by  its 
coordinates  x  =  (x1,x2,  x3)  and  velocity  components  c  =  (cy,  c2,  c3).  The  state  of  a  gas  is 
characterized  by  the  one-particle  distribution  function  /(x,  c,  t)  such  that  /(x,  c,  Qdxdc 
gives  at  time  t,  the  number  of  molecules  in  the  volume  element  with  position  vectors 
within  the  range  x  and  x  +  dx  and  with  velocity  vectors  within  the  range  c  and  c  +  dc. 

An  elastic  binary  collision  is  characterized  by  the  conservation  laws  of  momentum  and 
kinetic  energy,  namely, 

1111 

i  /  i  /  -l2i±2-l/2,-l/2  /i\ 

me  +  me i  =  me  +  mc1,  -me  2mCl  =  2mC  ^  2mCl  ’  (1) 

where  (c,  Ci)  refer  to  pre-collisional  velocities  and  (c^cQ  to  post-collisional  velocities  of 
two  molecules.  The  subindex  1  is  used  in  order  to  distinguish  two  identical  molecules  that 
participate  in  the  collision.  The  relative  velocity  is  denoted  by  g  =  (q  —  c  and  the  energy 
conservation  law  implies  that  g  =  g'. 

The  space-time  evolution  of  the  one-particle  distribution  function  in  the  phase  space 
is  governed  by  the  Boltzmann  equation,  which  in  the  absence  of  external  forces,  reads 

^  +  Ci<lr~  =  f  (f if' ~  fif)gbdbdedcl.  (2) 

We  note  that  it  is  a  non-linear  integro-differential  equation  for  /(x,  c,  t).  The  right-hand 
side  of  the  Boltzmann  equation  is  related  to  the  collisions  of  the  molecules  through  the 
product  of  two  distribution  functions.  The  relative  motion  of  the  molecules  is  characterized 
by  the  impact  parameter  0  <  b  <  oo  and  by  the  azimuthal  angle  0  <  £  <  27t  (see 
Figure  1).  Furthermore,  the  following  abbreviations  were  introduced  f  =  /(x,  c',t),  f[  = 
/(x,c'i,i),  /  =  /(x,c,t),  fi  =  /(x,d,Q. 

The  multiplication  of  the  Boltzmann  equation  (2)  by  an  arbitrary  function  -0  = 
0(x,  c,  t)  and  the  integration  of  the  resulting  equation  over  all  values  of  the  velocity 
components  c  leads  to 


d_ 

dt 


0/  dc  + 


d_ 

dxi 


iicj  dc  - 


90  90 

dt  Cl  dxi 


f  dc 


fw-mfsK***** c 


\ /(0i  +  0 -0'i -0/)(/i/'  -  fi^gbdbdedc^c,  (3) 


where  the  two  equalities  in  the  right-hand  side  of  the  above  equation  where  obtained  from 
the  symmetry  properties  of  the  collision  operator. 
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Figure  1:  Geometry  of  a  binary  collision. 


From  the  analysis  of  the  last  equality  in  the  right-hand  side  of  (3)  we  may  infer  that 
it  vanishes  for  any  kind  of  distribution  function  /  when  ipi  +  t/j  —  ^jj'x  +  A  function  -0 
which  fulfills  such  condition  is  called  a  summational  invariant.  A  summational  invariant 
is  given  by  a  linear  combination  of  the  conservative  quantities:  mass,  linear  momentum 
and  energy,  and  is  expressed  by  ^(c)  =  A  +  B  •  c  +  Dc 2  where  A  and  D  are  two  scalar 
functions  and  B  a  vectorial  function,  all  of  them  being  independent  of  c. 

In  kinetic  theory  a  macroscopic  state  of  a  gas  is  characterized  by  quantities  that  are 
defined  in  terms  of  the  distribution  function  /(x,  c,t).  Firstly,  based  on  the  microscopic 
quantities  of  the  gas  molecules  as  mass  m,  linear  momentum  mq  and  energy  me2  / 2  we 
define  the  mass  density  g,  the  momentum  density  gvt  and  the  total  energy  density  gu  of 
the  gas,  namely, 


g  —  mf(x.,c,t)dc,  gvt  —  mcif(x.,c,t)dc,  gu  = 


me 


/(x,  c,t)  dc.  (4) 


If  we  substitute  the  molecular  velocity  in  (4)3  by  the  peculiar  velocity  C*  =  q  —  vt  -  such 
that  f  mCif  dc  =  0  -  we  obtain 


-L  O 

gu  =  -gu  +  ge, 


where 


mC2/(x,  c,  t)  dc. 


(5) 


Hence,  the  total  energy  density  of  the  gas  is  given  by  a  sum  of  its  kinetic  energy  density 
gv2/ 2  and  its  internal  energy  density  ge. 

We  define  the  moment  of  the  distribution  function  of  order  N  by 


Phi2...iN  =  J  rnChCi2 . . .  CiNf(x,  c,  t)  dc ,  (6) 

which  represents  a  symmetric  tensor  of  order  N  with  (A^+l)(A^+2)/2  distinct  components. 

The  zeroth  moment  represents  the  mass  density,  the  first  moment  vanishes  and  the 
second  moment,  known  as  the  pressure  tensor,  reads 


Pij 


J  mCiCjf(x.,  c,  t)dc. 


(7) 
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We  denote  by  p^j)  the  pressure  deviator,  which  represents  the  traceless  part  of  the  pressure 
tensor.  It  is  given  by 


P{ij) 


Pij  ~ 


with 


md2/(x,  c,  t)  dc, 


(8) 


representing  the  hydrostatic  pressure  of  the  gas.  The  internal  energy  density  ge ,  given 
by  (5)2  is  related  to  the  hydrostatic  pressure  by  p  =  2ge/3.  The  equation  of  state  of  an 
ideal  gas  is  given  by  p  =  nkT  where  k ,  n  =  g/m  and  T  denote  the  Boltzmann  constant, 
the  particle  number  density  and  the  temperature  of  the  gas,  respectively.  Hence,  we  can 
obtain  the  following  expression  for  the  temperature  of  a  monatomic  gas  written  in  terms 
of  the  distribution  function: 


T  = 


P_ 

nk 


2  m  1 

3  k  3nk 


mC2/(x,  c,  t)  dc. 


(9) 


The  heat  flux  vector  is  defined  as  the  contracted  third  order  moment 


Qi  = 


1 

2  PiP 


J  mC2Cif(x,  c,  t)  dc, 


(10) 


and  the  moments  of  higher  order  do  not  have  specific  proper  names. 

We  obtain  the  balance  equations  for  the  moments  of  the  distribution  function  from 
the  transfer  equation  (3)  by  choosing  if{x,c ,t)  equal  to: 

(i)  Balance  of  mass  density  (if  =  m): 


dg  dgvj 
dt  dxi 


0, 


(11) 


(ii)  Balance  of  momentum  density  (if  =  mcf): 


dgvi  d 

-5r  +  g-(^+P«)  =  0, 


(12) 


(iii)  Balance  of  total  energy  density  (if  =  me2/ 2): 


d 

\  (  1  ,\1 

d 

X  (  1  I 

Wt 

K£+2”)J 

+  dxi 

+  2V  )Vi  +  qi+Piivi 

=  0, 


1112-IN  ,  d  ,  \  _  &PiN)k 

dt  +  dxk[Pt1t2...tNk+P*1V2...*ltVk)  ^/'hlO-CV  1  QXk 

dviN ) 


(iv)  Balance  of  Nth  order  moment  (if  =  rnC^ C(2 . . .  CiN): 
dp 

w.  -  “I"  r^VPi  1  io...i\rk  -t -Pi  -\io  ...iKT  Vk)  - 

Q 

f)lh..\ 

+  Npk( 


=  p 

(U12...JAT-1  Qx  1 


(13) 


(14) 


The  above  equation  was  obtained  by  eliminating  the  time  derivative  of  the  hydrodynamic 
velocity  Vi  by  the  use  of  (12).  The  parenthesis  around  the  indexes  indicate  a  sum  over  all 
AM  permutations  of  these  indexes  divided  by  N\.  Furthermore,  the  production  term  due 
to  the  molecular  collisions  reads 


P 

±  i\%2 ...in 


=  hn(C'ilC'i2...C'iN-CilCi2 


CiN)f  figbdb  de  dc  dc\ 


(15) 
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If  we  subtract  from  the  balance  equation  of  the  total  energy  density  (13)  the  mo¬ 
mentum  balance  equation  (12)  multiplied  by  the  bulk  velocity  Vi,  we  obtain  the  balance 
equation  for  the  internal  energy  density,  namely, 

f+±.t(e£Vi  +  q,)+Ptj^=0.  (16) 

3  Chapman-Enskog  Method 

First  let  us  introduce  the  dimensionless  variables 


t*  =  t/tc,  x*  =  x/Lc,  c*  =  c/c,  g*  =  g /g,  b *  =  b/ d,  £*  =  e/n,  f*  =  f/fc, 

(17) 

where  d  denote  the  molecular  diameter,  c  =  y/8 kT /mir  the  mean  thermal  velocity,  and  g  = 
\/2c  the  mean  relative  velocity.  Furthermore,  the  characteristic  value  for  the  distribution 
function  fc  is  equal  to  n/c3  where  n  represents  the  particle  number  density  of  the  gas  and 
tc  and  Lc  are  a  representative  time  and  length  that  characterize  the  fluid  flow.  By  taking 
into  account  the  dimensionless  variables  (17),  the  Boltzmann  equation  (2)  can  be  written 
as 

Srw +c4 § = 7  / (/r/"  ■ nr)g' v  db'  **  (is) 

In  the  above  equation  we  have  introduced  the  Strouhal  number  Sr  and  the  Knudsen 
number  Kn  which  are  dehned  by 


x*  =  x/Lc,  c*  =  c/c,  g*  =  g /g,  b *  =  6/d,  £*  =  e/7r, 


c  T'c 

Sr  =  — , 

Clr 


Kn  =  — , 


where  /  =  l/\/27rd 2n  is  the  mean  free  path. 

The  Knudsen  number  is  related  to  the  degree  of  rarefaction  of  a  gas.  When  Kn  <C  1 
the  molecular  collisions  are  very  important,  the  distribution  function  is  determined  by  the 
collision  term  of  the  Boltzmann  equation  and  the  gas  is  described  by  a  continuum  regime. 
Otherwise,  when  Kn  3>  1  the  molecular  collisions  become  negligible,  the  distribution 
function  is  determined  via  a  collisionless  Boltzmann  equation  and  the  regime  of  the  gas 
is  known  as  free  molecular  flow. 

In  the  method  of  Chapman-Enskog  the  distribution  function  /  is  expanded  as 

OO 

/  =  /(°)  +  A fW  +  A 2/(2)  +  • ' '  =  A7(r),  (20) 

r=0 

i.e.,  in  power  series  of  a  parameter  A  which  is  of  order  of  the  Knudsen  number.  The 
distribution  functions  f(°\  /  A)  and  represent  the  first,  second  and  third  approximation 
to  the  distribution  function,  and  so  on.  The  parameter  A  can  be  set  later  equal  to  unity, 
so  that  the  proper  dimensions  of  the  Boltzmann  equation  are  restored.  Moreover,  the 
approximations  must  satisfy  the  constraints 


ifjf^dc  =  0,  Vr  >  1, 
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where  i/>  represents  the  summational  invariants  m,  mci  (or  mCj)  and  me2/ 2  (or  mC2/ 2). 

If  we  insert  the  series  expansion  of  the  distribution  function  (20)  into  the  dehnitions 
of  the  pressure  tensor  and  of  the  heat  flux  vector  we  obtain 

-j  n  n  OO  OO 

Pij  =  Sij  -  /  mC2f^dc  +  /  mC^Cj)  ]T  Arf^dc  =  P6tj  +  ^  Arp{g ,  (22) 

J  r= 1  r=  1 

OO  OO 

jCfiC,  J2  A7<r)*  =  A'#>.  (23) 

r=l  r= 1 


The  parameter  A  is  introduced  into  the  collision  term  of  the  Boltzmann  equation  (2), 
yielding 

Vf+c€riQifJ)’  (24) 

where  V  =  d/dt  +  Vid/dxi  denotes  the  material  time  derivative.  In  (24)  the  collision  term 
was  expressed  in  terms  of  an  integral  for  a  bilinear  quantity 

Q(F,G)  =  ^  I (F[G'  +  F'G'1-F1G-FG1)gbdbd£dc1.  (25) 

The  material  time  derivative,  like  the  distribution  function,  is  also  expanded  in  power 
series  as 

OO 

V  =  V0  +  AV1  +  A2V2  +  ...  =  J2  Ar£V,  (26) 

r= 1 

and  the  insertion  of  the  expansions  (22),  (23)  and  (26)  into  (11),  (12)  and  (16)  lead  to 
the  following  decomposition  of  the  balance  equations 


d  Vi 

V0p  +  — 

OXi 

=  0, 

-i 

ro 

O 

'< 

IV 

(27) 

^  dp 

QF>o  Vi  +  — 
OXi 

=  0, 

dp(r). 

gVrVi+  {j)  =0  (Vr>  1), 

dxj 

(28) 

3  d  Vi 

-nkVrT  +  ^  +4,)&,=°  (Vr£1)' 

(29) 

=  0, 

We  can  now  obtain  the  integral  equations  for  the  approximations  of  the  distribution 
function  by  the  inserting  the  expansions  (20)  and  (26)  into  the  Boltzmann  equation  (24) 
and  equating  equal  powers  of  A.  Hence,  it  follows 

Q(/(0),/(0))  =  0,  (30) 

2  Q(/(0),/(1))  =  Po/(0)  +  a^,  (31) 

2  Q(/(0),/(2))  +  Q(/(1),/(1))  =  Vof^+VJ^+Q^,  (32) 

and  so  on.  The  above  equations  represent  the  three  first  integral  equations  for  /(°), 
and  f(2\  respectively. 
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For  the  determination  of  the  Erst  approximation  to  the  distribution  function  /(0^  from 
the  integral  equation  (30)  we  note  that  its  solution  is  given  by  y'(o)y4°)  =  /(°)/|°)  or 
In  /'(°)  +  ln/j(°^  =  In  +  In  /,'0) .  Hence,  In  is  a  summation  invariant  so  that  it  must 
be  represented  by  In  /(0)  =  A  +  B  •  C  +  DC 2 .  If  we  insert  this  representation  for  / into 
the  definitions  of  the  mass  density  (4)!,  hydrodynamic  velocity  (4)2  and  temperature  (9) 
we  obtain  the  Maxwellian  distribution  function,  namely, 

3 

/(0)  =  n  (-Y  e-^2  where  (33) 

\7 t)  2 kT 

For  the  determination  of  the  second  approximation  / ^  we  write  =  f^(f>  and  the 
integral  equation  (31)  reduces  to 

23(/<0|,/(0V)  =  /<“,{l(^-*)C,g  +  2/3C,C,g} 

=  1[4>]  =  I  /(0)/i(0)  (0i  +  4>'  ~4>i  ~4>)gbdbde  dcq,  (34) 


by  eliminating  the  material  time  derivatives  through  the  use  of  (27)  i,  (28)  i  and  (29)  i. 

The  general  solution  of  the  integral  equation  (34)  is  given  by  a  sum  of  the  solution 
of  the  homogeneous  integral  equation  plus  a  particular  solution  of  it,  i.e.,  (f)  =  4>h  +  4>P- 
The  solution  of  the  homogeneous  integral  equation  I[(f>h]  =  0  is  a  summational  invariant 
c i>h  =  ai  +  a^Cj.  +  aaC2,  where  aq  and  a3  are  scalar  functions,  while  o;2  is  a  vector  function, 
all  of  them  do  not  depend  on  the  peculiar  velocity  C. 

The  particular  solution  of  the  integral  equation  is  expressed  as  a  linear  combination 
of  the  temperature  gradient  and  of  the  velocity  gradient  deviator,  namely, 


Ai  DT 
T  dxi 


-2/TR 


dv(j 

dxj) 


(35) 


where  Ai  and  BtJ  denote  a  vector  function  and  a  tensor  function  of  (C,  g,  T ),  respectively. 
It  is  easy  to  show  that  the  representations  of  the  vector  and  tensor  functions  in  terms 
of  the  peculiar  velocity  C  are  given  by  Ai  =  A*Ci  and  Bij  =  BCiCj,  where  the  scalar 
functions  A*  and  B  depend  only  on  (C2,  g,  T ). 

Hence,  the  solution  of  the  non-homogeneous  integral  equation  (34)  is  given  by 

0  =  -A’  ^  ^  -  2pBCiCj^  +  oq  +  a2rCr  +  as^2.  (36) 

The  second  approximation  fA  =  must  satisfy  the  constraints  (21)  which  imply 
the  following  restrictions 


3  kT 

O' i  H - 0(3  —  U, 

m 


5  kT 

o i  H - o3  —  U, 

m 


A*  d T 

~Ydxj 


/<0>dc  =  0. 


(37) 


From  the  two  first  equations  we  may  infer  that  aq  =  o3  =  0,  whereas  the  last  equation 
implies  that  o2  must  be  proportional  to  dT/dxj.  By  writing  the  proportionality  factor  of 
o2  as  o/T,  the  constraint  (37)3  reduces  to 


(38) 
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with  A  =  A*  —  a  being  a  new  scalar  function  of  ( C 2,  g,  T ). 

Hence,  from  the  above  results  we  may  write  the  distribution  function  (20)  -  up  to  the 
second  approximation  and  with  A  =  1  -  as 

/  =  /<“»(!  +  0)  =  /(0)  {l  -  ^C~  -  2/3BCiC,|'|}  .  (39) 

Now  the  scalar  coefficients  A  and  B  in  the  distribution  function  (39)  can  be  obtained 
as  solutions  of  the  integral  equations  that  follows  from  (34),  namely, 


/(0) 


C,  =  -I\AC,],  fmC(,Cj}  =  -I[  BCjiCj,]. 


(40) 


In  order  to  solve  the  integral  equations  (40)  the  scalar  coefficients  A(C2,  q,  T)  and 
B{C2,q,T)  are  expanded  in  series  of  Sonine  polynomials.  The  Sonine  polynomials  of 
order  n  and  index  (/  +  1/2)  in  the  variable  f3C2  =  mC2 /2kT  are  dehned  by 


siiUec2)  =  J2 


T(n  +  l  +  3/2) 


(-/ K2)k , 


k= 0 


k\(n  —  k)\T(k  +  l  +  3/2) 


(41) 


and  obey  the  orthogonality  conditions 

/»oo 

/  e-^C™s£\rSPC2)S<l\/2(l3C 


2)P,+3,2dC  =  T("+.(,+  3/2b„p.  (42) 


n\ 


where  T(rt)  denotes  the  gamma  function.  The  Erst  three  Sonine  polynomials  read 


Si%(P°2)  =  1,  S«/2CSC2)  =  l  +  pc2,  (43) 

s£/2(/3C2)  =  l(‘+l){‘+l)-(‘+l)PC2  +  lp2C2.  (44) 

The  integral  equations  (40)  in  terms  of  the  Sonine  polynomials  become 

f«»4%(pC1)Ci=X[ACi],  (PC2)C(lCj}  =  -X  [BCtfCfi]  .  (45) 

Moreover,  the  expansions  of  the  scalar  coefficients  M/C2,  g,  T )  and  B(C2,  g,  T )  in  series  of 
Sonine  polynomials  are  written  as 


A(C2,e,T)  =  -Y/a^S^(PC2),  B(C2,p,T)  =  J^b^S^PC2).  (46) 

r=0  r=0 


In  the  above  equations,  the  scalar  coefficients  a ^  and  b ^  are  only  functions  of  (g,T), 
and  from  this  point  on,  the  dependence  of  the  Sonine  polynomials  in  the  variable  (3C 2  is 
omitted. 

The  coefficient  A  must  satisfy  the  constraint  (38),  hence  by  using  its  representation 
(46)  i,  we  get 

/°°  /  q  \  §  00  r°°  o 

C2f»  J2  a^S^dc  =  4™  “W  /  S^S^e-^dC  =  - -<P\ 

r=0  b  '  r=0  ^ 

(47) 
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clue  to  (42)  and  (43) i.  From  (47)  it  follows  that  =  0  and  the  scalar  coefficient  A 


reduces  to 

OO 

A(C2,g,T)  =  -J2^r)Sir/)2.  (48) 

r= 1 

Now  the  insertion  of  the  expressions  (48)  and  (46)2  into  their  respective  integral  equations 
(45),  yields 


£■ 

r= 1 


i{r)l 


q(r)  f  i 
°3/2 


b(r)l 


r= 0 


(49) 


In  order  to  determine  the  coefficients  a,A  e  bA  from  (49)  we  proceeds  as  follows.  First 
we  multiply  the  integral  equation  (49)  i  by  (3Sy2Ci  and  integrate  the  resulting  equation 
over  all  values  of  C.  By  using  (42)  we  obtain 


15  n 


<$(1’r)  = 


S=1 


a(r’*Vs), 


where 


(50) 


In  the  above  equation  represents  Kronecker’s  symbol.  Following  the  same  method¬ 
ology,  we  multiply  the  integral  equation  (49)2  by  fAS^C^Cj) ,  integrate  the  resulting 
equation  and  get 


--^(0’r)  =  J2P{r,s)b{s\  where  /3(r’s) 

'  s=0 


S^C^HS^C^dC.  (51) 


The  integrals  in  (50)2  and  (51)2  depend  on  the  molecular  interaction  potential  and  on  the 
fields  ( g,T ).  Moreover,  from  the  condition  f  (fl[(p\dc  =  f  (j)X[(p\dc  -  which  is  valid  for  all 
arbitrary  functions  0(x,  c,  t)  and  <^(x,  c,  t)  -  we  may  infer  that 

O.AA  =  aAA,  p(r,s)  =  p(s,r)  and  a(r,0)  =  a(0,r)  =  Q.  (52) 


We  conclude  that  (50)i  and  (51)i  represent  an  infinite  system  of  algebraic  equations 
for  the  coefficients  and  An\  respectively. 

To  obtain  the  constitutive  equations  for  the  pressure  tensor  and  for  the  heat  flux 
vector  we  proceed  as  follows.  First,  we  insert  the  non-equilibrium  distribution  function 
(39)  together  with  the  representations  of  A  and  B  -  given  by  (48)  and  (46)2,  respectively 
-  into  the  definitions  of  the  pressure  tensor  and  heat  flux  vector  and  get 


dvr. 


Pij  =  I  mCiCjfdc  =  p  Sij  -  I  2 mpCiCjS^fW  ^  S{5%C{pCq)^dc,  (53) 


r= 0 


m 


Qi  =  I  y C2CJdc  = 


m 


2  Tf3 


-/*  -  411  C‘f(0,Za<r,sf/AA:d c-  <54> 


In  (54)  we  have  used  the  relationship  f3C2  = 


5  c(°)  _  o' 

2°3/2  J3/2 


dT 

r=l  3XP 

(1)' 


which  follows  from  (43)  2. 


Next,  the  substitution  of  the  Maxwellian  distribution  function  (33)  into  (53)  and  (54)  and 


RTO-EN-AVT-194 


9-11 


3  CHAPMAN-EN SKOG  METHOD 


the  subsequent  integration  of  the  resulting  equations  over  all  values  of  C,  yield 
Pij  =  pSij  -  2 n~ 


dvp 
ldxj} ’ 

where 

u  =  — 6(0) 

1  2(3  ’ 

(55) 

dT 

A  — - , 

where 

A  = 

(56) 

dxi 

4m  (3 

thanks  to  the  orthogonality  condition  of  the  Sonine  polynomials  (42).  Equations  (55) i 
and  (56) i  are  the  mathematical  expressions  for  the  laws  of  Navier-Stokes  and  Fourier, 
with  (i  and  A  representing  the  coefficients  of  shear  viscosity  and  thermal  conductivity, 
respectively. 

We  note  from  (55)2  and  (56)2  that  the  transport  coefficients  of  shear  viscosity  and 
thermal  conductivity  are  functions  only  of  the  coefficients  a/ '1  and  b^\  respectively.  We 
can  obtain  these  coefficients  from  the  system  of  equations  (50)i  and  (51)i,  which  can  be 
written  as 


aw  =  lim  kk. 

n^oo  Ann 


&(°)  =  lim 

n-HX>  Bn  n 


(57) 


where 


A'nn  =  det 

.  0  ^ 

c.41’2)  .. 

A2^  .. 

.  a(bn) 

.  A2A 

)  d.6t 

«(bi)  .. 

.  cdbA 

0 

a(na)  . . 

(58) 


B'nn  =  det 

5  n 

2  p2 

0 

pw  .. 

/3(bd  .. 

.  p<-  bn) 

7  $nn  det 

.  p^A 

0 

/3(n’b  . . 

Q(n,n) 

/5(n,°)  _  _ 

j3{n,n) 

(59) 


The  coefficients  a-b  and  b ^  follow  from  (58)  and  (59)  through  a  method  of  successive 
approximations,  where  the  approximation  of  order  p  for  the  coefficients  a d)  and  b - 
denoted  by  [ad)]p  and  [6®]p  -  is  given  by  the  sub-determinant  of  order  p.  The  first  two 
approximations  for  the  coefficients  a- b  e  b ^  are 


[«(1)]i  = 


15  n 


a(1)]2  = 


15  n 


4  /32  cdbi)  ’ 

a(2.2) 


=  I"  1 


2  /32  pm  ’ 


2,  r']2=~ 


^(bb 


4  P2  a(i-i)a(2,2)  _  a(i,2)2  ’  L“  JZ  2  (32  pmp^A)  _  ^(o,i)2 ' 
The  integrals  cobW  and  pr,s )  are  given  in  terms  of  the  collision  integrals 


n{l’r)  = 


POO  POO 


'0  JO 


e  '>,2y2r+3  (l  —  cosz  x)  b  db  dp, 


(60) 

(61) 

(62) 


where  \  =  arccos(g'  •  g / g2)  is  the  scattering  angle  and  7  =  \J  (3/2g  represents  a  dimen¬ 
sionless  relative  velocity. 

The  determination  of  the  integrals  a d’s)  and  /?(r>s)  is  not  an  easy  task.  As  an  example, 
let  us  compute  the  integral  /^0,°).  For  its  calculation  a  change  of  the  integration  variables 
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is  performed  by  replacing  the  velocities  C  and  Ci  by  the  relative  velocity  g  =  Ci  —  C  and 
by  the  center  of  mass  velocity  G  =  (Ci  +  C)/2.  The  linear  momentum  conservation  law 
implies  that  G'  =  G  and  the  Jacobian  of  the  transformation  of  the  volume  elements  has 
modulus  equal  to  one  so  that  dCdCi  =  dgd.G.  In  terms  of  the  new  variables  the  integral 
reads 

=  /  C(iCy)I[C(iC3)]<ic  =  d  (ty  j e-WO'e-S'' 
x  (g>,Gj)  -  Gi.tj,  +  (7 ,.,/,)  (jjjj',  -  7,'7;  )  gbdldedgdG.  (63) 
The  integration  of  (63)  with  respect  to  G  and  £  leads  to 

3 

/5(0,o)  =  J  e~^2  (cos2  x  -  1  )g5bdb  dg,  (64) 

If  we  introduce  the  dimensionless  relative  velocity  7  =  y//3/2g,  the  integration  of  the 
above  equation  with  respect  to  the  angles  of  g  leads  to 

/?<”•»>  =  -V?n2  (H  "  Ji(2’2),  (65) 

where  denotes  the  collision  integral  defined  by  (62). 

By  introducing  the  integrals  air,s)  =  J^at^r’s\  and  =  n2^ftr,s\  we  can  obtain 
in  the  same  way  the  following  expressions  in  terms  of  the  collision  integrals  QSl,r) 

a<M)=4SiM  o<u>  =  7«(2'2>  -  2f!<2X  aM  =  ^!JM_7SiM  +  S1(2,4)j  (66) 
=  7U),  7“4)  =  ol1'21.  7U)  =  ^!2(2’2>-7!J<2’3>+!J‘2’4>  = 

(67) 

As  was  pointed  out  the  coefficients  b 1°)  and  a-d  are  determined  through  a  method  of 
successive  approximations,  so  are  the  transport  coefficients  /i  and  A.  Hence,  it  follows 
from  (60),  (61),  (55)2  and  (56)2  that  the  two  first  approximations  to  the  coefficients  of 
shear  viscosity  and  thermal  conductivity  read 


Mi  =  j 


mkT  1 


7T 


(0,0)  ’ 


mkT 


vr  ^,0^1,1) 


M1,2)2’ 


[A]2  = 


75  k 
16m 


[A]i  = 

mkT 
7 r 


75  k  mkT  1 


16  m 


7 r 


(2,2) 

tt* 


a 


(i,i)’ 


(1,1)  (2,2) 
a*  oil 


(1,2) 

a* 


2  • 


(68) 

(69) 


Now  from  (68)  and  (69)  we  can  build  the  ratios 

Mi  5  [A],  /  tt72>2  /  <r)a  y1  M2 

Mi  2”’  [A],  (  )  {  I3lm  tf’1) )  Mi' 


(70) 
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We  may  infer  from  (70) i  that  the  ratio  of  the  Erst  approximations  to  the  coefficients  of 
thermal  conductivity  and  shear  viscosity  of  monatomic  ideal  gases  is  equal  to  5cv/2  for 
all  types  of  spherically  symmetrical  molecular  interaction  potentials,  where  cv  =  3k /2m  is 
the  specific  heat  at  constant  volume.  Furthermore,  from  the  relationship  (70)2  it  follows 
that  [A]2/[/z]2  >  [A] i/[/i] i  =  5c„/2,  i.e.,  the  ratio  of  the  second  approximations  is  larger 
than  the  ratio  of  the  first  ones. 

For  a  hard-sphere  potential  the  impact  parameter  is  given  by  b  =  dcos(y/2)  with 
0  <  x  <  7r,  so  that  the  collision  integral  rh2,ri  becomes 


roo  i 

=  /  e^yr+3-d2d7 
Jo  3 

Hence,  it  follows  from  (66)  and  (67)  that 
aa,i)  =  =  4^  aM  =  „(»,.)  =  _d*, 


(71) 


205 

~Y2 


d2, (72) 


and  the  two  first  approximations  (68)  and  (69)  to  the  coefficients  of  shear  viscosity  and 
thermal  conductivity  for  the  hard-sphere  potential  read 


r  .  5  ImkT  [Ah  5  [uL 

■  \i=2c- 


205 

202 


1.014  851, 


[Ak 

[A], 


45 

44 


1.022  727.(73) 


We  note  that  the  second  approximations  to  the  transport  coefficients  is  a  small  correction 
to  the  first  ones,  so  is  the  ratio  [A]2/[/i]2  ~  1.007 761[A]i/[/x]i. 

The  fifth  approximations  to  the  transport  coefficients  for  hard-sphere  potential  are 


j^A  ~  1.016  027,  «  1.025  197. 

Ni  [A]i 


4  Application  of  Chapman-Enskog  Method  to  Gran¬ 
ular  Gases 

The  mechanical  energy  of  a  gas  is  conserved  when  its  molecules  undergo  clastic  colli¬ 
sions  and,  in  this  case,  the  gas  relaxes  towards  an  equilibrium  state  characterized  by  a 
Maxwellian  distribution  function.  However,  the  inelastic  collisions  of  the  gas  molecules 
transform  the  translational  kinetic  energy  into  heat  and  the  mechanical  energy  lost  im¬ 
plies  a  temperature  decay  of  the  gas.  In  the  literature  gases  whose  molecules  undergo 
inelastic  collisions  are  known  as  granular  gases.  The  main  premisses  of  the  kinetic  theory 
of  granular  gases  are:  (a)  only  binary  collisions  of  hard  spherical  molecules  are  taken  into 
account  and  (b)  the  energy  lost  from  inelastic  collisions  is  small. 

Let  us  consider  the  encounter  of  two  identical  molecules  of  mass  m  diameter  d  pre- 
collisional  velocities  (c,cq)  and  post-collisional  velocities  (c',^).  The  momentum  conser¬ 
vation  law  mc+mci  =  mc'+mc)  holds  for  inelastic  collisions,  but  the  inelastic  encounters 
are  characterized  by  the  relationship  (g'  •  k)  =  —  a(g  •  k)  which  relates  the  pre-  and  post- 
collisional  relative  velocities  at  collision.  The  parameter  0  <  a  <  1  -  here  considered  as 
a  constant  -  refers  to  the  normal  restitution  coefficient  and  k  is  the  unit  vector  directed 
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along  the  line  which  joins  the  molecules  centers  and  pointing  from  center  of  the  molecule 
labeled  by  the  index  1  to  the  center  of  the  molecule  without  label.  The  component  of  the 
velocity  perpendicular  to  the  collision  vector  k  does  not  change  in  inelastic  collisions,  so 
that  k  x  g'  =  k  x  g.  From  the  momentum  conservation  law  it  follows  the  relationships 
which  give  the  post-collisional  velocities  in  terms  of  the  pre-collisional  ones,  as  well  as 
relationships  which  connect  the  relative  velocities  and  their  modulus: 

c'  =  c  +  k)k,  ci  =  C!  -  iy^(g-k)k,  (74) 

g/  =  g-(l  +  a)(g-k)k,  g'2  =  g2  -  (1  -a2)(g  •  k)2.  (75) 


In  the  case  of  elastic  collisions,  a  =  1  and  it  follows  the  conservation  of  the  kinetic  energy 
9  =  9'- 

A  restitution  collision  with  pre-collisional  velocities  (c*,c*)  that  corresponds  to  the 
post-collisional  velocities  (c,  Ci)  are  characterized  by  the  equations 

c*  =  c+  ^^(g  •  k)k,  c*  =  ci  -  t^(g-k)k,  (76) 

where  k*  =  — k  and  (g  •  k)  =  —  a(g*  •  k). 

For  the  determination  of  the  Boltzmann  equation  we  have  to  know  the  transformation 
of  the  volume  elements  dc\  dc*  =  \J\dc\dc,  where  the  modulus  of  the  Jacobian  is  given 
by  |  J\  =  1/a.  Hence,  it  follows  the  relationship 

(g*  •  k*)  dc*  dc\  =  -?(g  ’  k)  dcdci,  (77) 

cr 

and  the  Boltzmann  equation  for  granular  gases  reads 

%+*§/  =  I  (A/i*/'-/i/)  d2(g.k)dkrfCl  =  Q, (/,/).  (78) 

Above  we  have  introduced  the  bilinear  form 

Qj(F,G)=l-  J  (^F*G*  +  ^r2F*Gl-  F1G  -  FG^  d2  (g  •  k)  dkdc,.  (79) 


The  balance  equations  for  the  mass  density  g,  for  the  momentum  density  ov%  and 
for  the  specihc  internal  energy  £  =  3kT/2m  are  obtained  from  the  Boltzmann  equation 
(78)  by  multiplying  it  successively  by  m,  mci  and  mC2 / 2  and  integrating  the  resulting 
equations.  Hence,  it  follows 


Vg  +  g 


dvj 

dxi 


0, 


Vvi  + 


dpij 

dxj 


0, 


VT  + 


2 

■ink 


+  Pij 


dvj\ 

dxj) 


+  T( 


0.  (80) 


We  observe  that  the  balance  equation  for  the  temperature  (80)3  has  the  additional  term 
TQ  which  represents  the  energy  loss  due  to  inelastic  collisions.  The  coefficient  £,  known 
as  cooling  rate  of  the  granular  gas,  is  given  by 


d2m(l  —  a2) 
12  nkT 


/i/(g  •  k)3dk  dci  dc. 


(81) 
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It  is  easy  to  verify  from  the  above  equation  that  the  cooling  rate  vanishes  for  elastic 
collisions  where  a  —  1. 

The  methodology  used  to  determine  the  distribution  function  by  using  the  Chapman- 
Enskog  method  for  granular  gases  is  different  from  the  one  applied  to  a  monatomic  gas 
of  clastic  spherical  molecules,  since  we  have  to  take  into  account  that  there  exists  no 
equilibrium  state  characterized  by  a  Maxwellian  distribution  function.  The  first  difference, 
refers  to  the  decomposition  of  the  balance  equations,  since  they  are  written  as 


(J  U  ' 

V0n  =  0,  T>in  =  —n  — — (82) 

axi 

1  dv 

V0Vi  =  0,  V1  Vi  =  —  — .  (83) 

Q  OXi 

O T  rhi- 

VqT  =  -TC(0),  D1T  =  -TC(1) - —A  (84) 

3  OXi 

The  first  and  the  second  approximation  to  the  cooling  rate  in  (84)  read 

C(0’  =  d~’^2>  / 4°'/(0)(g  •  kfdkdc,  dc,  (85) 

c(1)  =  d2t-2)  J  (/!1)/(0)  +  /i(0>/(1))  (g  ■  k)Ukdci  dc.  (86) 


In  terms  of  the  parameter  A  we  write  the  Boltzmann  equation  for  granular  gases  (78) 


as 


1 

A 


Vf  +  Q 


df_ 

dxi 


{Qi  (/./). 


(87) 


instead  of  the  representation  (24),  indicating  that  the  material  time  derivative  and  the 
collision  term  are  of  same  order,  while  the  spatial  gradients  are  of  higher  order. 

The  two  first  integral  equations  for  the  determination  of  / and  are  obtained  by 
inserting  the  expansions  (20)  and  (26)  into  the  Boltzmann  equation  (87)  and  equating  the 
equal  powers  of  A,  yielding 


Vo /(0)  =  Q/(/(0),/(0)),  VJ^  +  Vo /(1)  +  =  2  Q/(/(1),/(0)).  (88) 

In  order  to  determine  the  first  approximation  of  the  distribution  function  the  integral 
equation  (88)  i  is  written  as 


-n  (0|^71  =  Q/(/<0),/,0)), 


(89) 


through  the  elimination  of  the  material  time  derivatives  T> o  by  using  (82)  i,  (83)  i  and 
(84)^  We  note  that  the  solution  of  the  integral  equation  (89)  for  the  distribution  function 
/D)  is  not  a  Maxwellian,  but  we  may  write  as  a  Maxwellian  plus  an  expansion  in 
Sonine  polynomials  of  the  peculiar  velocity,  i.e., 


/(0)  =  !m 


1  +  5>„S<1”)(/3C2) 

'  2 

71—  1 


where  —  e  /3c2,  (90) 
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denotes  the  Maxwellian  distribution  function.  The  scalar  coefficients  an  do  not  depend 
on  the  peculiar  velocity  C  and  we  suppose  that  this  series  expansion  converges  for  large 
values  of  the  peculiar  velocity.  We  may  approximate  the  representation  of  the  distribution 
function  (90)  as 


/(0) 

S{'\(3C2) 


fM  l  +  aiS{l\pC2)  +  a2Sy?\pC2) 


(2)/ 


where 


S?\l3C 2)  =  0  -  \/)C2  +  l-fC\ 

2  O  A  A 


(91) 

(92) 


If  we  insert  the  representation  (91)  into  the  definition  of  the  temperature  (9)  and  integrate 
the  resulting  equation  we  get  that  (i\  =  0. 

For  the  determination  of  the  coefficient  a2  we  proceed  as  follows.  First  the  product  of 
the  distribution  functions  is  written  as 


1  +  02 


15 

T 


+ 


P 

2 


^2G4  +  G2g2  + 


+  2(G  •  g) 


(93) 


by  neglecting  all  product  of  the  coefficient  a2  and  changing  the  velocity  variables  (C,  Cx) 
by  the  relative  and  the  center  of  mass  velocities  (g,  G).  Next,  we  multiply  the  integral 
equation  (89)  by  an  arbitrary  function  of  the  peculiar  velocity  il>(C2)  and  integrate  the 
resulting  equation  over  all  values  of  c,  yielding 

-  J  ?Hc2)rc<0>+Pdc  =  \  J  p(c? )  +  v(c°) 

-  l/.(C12)-V(C'2)]/1<“,/<",d2(g-k)<lk<lc1<lc.  (94) 

If  we  choose  C 2)  =  1  and  ^(C2)  =  C 2  in  (94)  the  integration  of  the  resulting  equations 

lead  to  identities.  However,  if  we  choose  ■?/?( C 2)  =  C4  in  (94)  and  use  the  relationship 

Cf  +  C'4  -  Cf  -Cl  =  2(a  +  l)2(g  ■  k)2(G  •  k)2  +  L+iL(g  .  k)4 


+(a2  -  l)(g  ■  k)2G2  +  PUilfg  .  k)292  -  4 (a  +  l)(g  •  k)(G  •  k)(G  ■  g),  (95) 

we  obtain  through  the  integration  of  the  resulting  equation  that  the  coefficient  <22  has  the 
form 


16(1  — a)(l  — 2ct2) 

81  —  17a  +  30o;2(l  —  a)  ’ 


(96) 


We  may  infer  from  the  above  equation  that  the  coefficient  <22  vanishes  for  elastic  collisions 
where  a  —  1. 

From  the  knowledge  of  the  coefficient  a2,  the  distribution  function  (91)  is  written  as 


/(0)  =  Sm 


16(1  —  q)(l  —  2q2)  (2)  2 

81  -  17a  +  30a2(l  -  a)  5  U 


(97) 
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Figure  2:  Left  frame:  dimensionless  distribution  function  i: /n/33^2  as  a  function  of 
x  =  PC2  for  a  =  1  and  0.8.  Right  frame:  cooling  rate  C^/nd2  versus  the  temperature 
and  the  normal  restitution  coefficient. 


If  we  substitute  the  distribution  function  (97)  into  the  expression  for  the  cooling  rate 
(85)  and  integrate  the  resulting  equation  we  get  that 


c<°>=Ud2(i 

3  V  m 


1  + 


3(l-a)(l-2a2) 

81  —  17a  +  30a2(l  —  a)_ 


(98) 


In  the  left  frame  of  Figure  2  it  is  represented  the  dimensionless  distribution  function 
7r3/2y(o)^n^3/2  ag  a  function  0f  x  =  /3C2  and  a.  We  may  infer  from  this  figure  that  the 
peak  of  the  curve  decreases  by  decreasing  the  normal  restitution  coefficient.  In  the  right 
frame  of  the  Figure  2  it  is  plotted  the  cooling  rate  £(0)/nd2  versus  the  temperature  and 
the  normal  restitution  coefficient.  We  note  that  the  cooling  rate  increases  by  decreasing 
the  normal  restitution  coefficient  and  by  increasing  the  temperature.  Furthermore,  the 
cooling  rate  is  zero  when  a  =  1,  i.e. ,  for  clastic  collisions  of  the  molecules. 

The  determination  of  the  second  approximation  to  the  distribution  function  from 
the  integral  equation  (88)2  is  more  involved  and  the  details  of  such  kind  of  calculation  will 
not  be  given  here.  From  this  integral  equation  we  may  conclude  that  is  a  function 
of  the  thermodynamic  forces:  deviator  of  the  velocity  gradient  dvu/dxj \,  gradient  of  the 
particle  number  density  dn/dxi  and  gradient  of  temperature  dT/dxi.  Hence,  we  may 
write  the  second  approximation  as 


/(1)  =  /m 


7i^1)(C'2)G4 


<91nT 

dxi 


+  72  +  73  S)?(C2)Ct 


dvt 


(!)/ 


dx 


3) 


d  In  n 
dxi 


(99) 


where  71,  72  and  73  are  scalar  coefficients  that  depend  on  n,  T  and  a.  If  we  solve  the 
integral  equation  (88)2  we  get  that  the  coefficients  of  the  distribution  function  are 
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given  by 


7i 


7s 


15  /  m  —15 

d2n(9  +  7a)(a  +  l)  V  Trfcf’  72  “  2d2n(13  -  a) (a  +  1) 
300(1  —  a)  /  m 

d2n(9  +  7a)(a  +  1)(19  -  3a)  V  ' 


(101) 


We  obtain  the  constitutive  equations  for  the  pressure  tensor  and  heat  flux  vector 
through  the  substitution  of  the  distribution  function  (99)  into  their  definitions  and  by 
integrating  the  resulting  equations.  Hence,  it  follows  the  laws  of  Navier- Stokes  and  Fourier 


Pij 


Qi 


J  mCiCjfdc  - 

J  jC2ajdc 


p  Sij  -  2/i 

dT 

dxi 


dvr 


dx 


j) 


■0T 


d  In  n 

dxi 


(102) 

(103) 


Above,  the  coefficients  of  shear  viscosity  // ,  thermal  conductivity  A  and  the  one  associated 
with  the  gradient  of  particle  number  density  i)  read 


d 

■0 


15  ImkT  75  k  ImkT 

2d2(13  —  a) (a  +  1)  V  tt  ’  “  2d2(9  +  7a)(a  +  1)  m  V  vr 

750(1  -a)  k  ImkT 

d2n(9  +  7a)  (a  +  1)(19  -  3a)  mV  vr 


(104) 

(105) 


We  note  that  the  coefficients  are  proportional  to  \/T  so  that  they  increase  by  increasing  of 
the  temperature.  With  respect  to  the  normal  restitution  coefficient  a,  we  may  infer  that 
the  coefficients  also  increase  by  decreasing  the  values  of  a.  Furthermore,  by  considering 
elastic  collisions  where  a  =  1  the  coefficients  of  shear  viscosity  and  thermal  conductivity 
reduce  to  the  Erst  approximations  (73)  1,2,  while  the  coefficient  i)  vanishes. 


5  Grad’s  Method  of  Moments 


Let  us  consider  a  thirteen  scalar  held  description  of  a  rarefied  gas  within  Grad’s  moment 
method.  In  this  case  the  fields  are  the  mass  density  g(x.,t),  the  hydrodynamic  velocity 
Vi(x,t),  the  pressure  tensor  ptJ(x,  t)  and  the  heat  flux  vector  qi(x.,t )  and  their  balance 
equations  read 


dg 

dt 

dgiq 

dt 

dPij 

dt 

dcp 

dt 


+ 

+ 

+ 

+ 


dgiq. 


=  0, 


dxi 

d  (gVjVj  +  Pij ) 

dxj 

d  0 Pijk  +  PijVk )  dvj 


=  0, 


dxk 

d  ( Qij  +  qjVj) 
dXi 


+  Pki 


dxi 


+  Pkj 


dvi 


=  Pi 


ij, 


+  Pijk 


dvi 


dxk  ^  dx 


dxk 

dvi  pki  dpkj 


g  dXj 


1  Prr  dpij 

2  g  dxj 


Qi- 


(106) 

(107) 

(108) 
(109) 
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We  note  that  the  two  first  equations  above  represent  the  balance  equations  of  mass  density 
(11)  and  moment  density  (12),  while  the  balance  equation  for  the  pressure  tensor  (108) 
was  obtained  from  (14)  by  taking  N  =  2.  The  balance  equation  for  the  heat  flux  vector 
(109)  follows  from  (14)  by  taking  N  =  3,  contracting  two  indices  and  dividing  the  resulting 
equation  by  2.  Moreover,  we  have  introduced  the  following  definitions  in  (108)  and  (109): 

Piik  =  f  mCiCjCkfdc,  Pij  =  Jm  (C'C'  -  ()([,)  ffglY.  (110) 

/TTJ  I  711 

-tfCiCjfdc,  Qi=Jj  (C'2C'  -  C2°i)  ff irfr-  (m) 

The  tensors  pl]k  and  qV}  represent  non-convective  fluxes  for  the  pressure  tensor  and  heat 
flux  vector,  respectively,  while  PtJ  and  Qi  refer  to  their  production  terms.  In  the  above 
equations  we  have  introduced  the  abbreviation  dT  =  gbdbdedcidc. 

If  we  analyze  the  system  of  partial  differential  equations  (106)  -  (109)  we  may  infer 
that  it  cannot  be  considered  a  system  of  held  equations  for  the  basic  fields  g,  vtl  pl3  and 
cp.  In  fact,  to  obtain  a  closed  system  of  differential  equations,  the  non-convective  fluxes 
Pijk,  Qij  and  the  production  terms  Pt] ,  Qi  must  be  expressed  in  terms  of  the  basic  fields, 
and  for  this  end  we  must  express  the  distribution  function  in  terms  of  the  thirteen  scalar 
fields  g,  vt,  pl3  and  qi. 

In  Grad’s  moment  method  the  distribution  function  is  expanded  in  series  of  tensorial 
Hermite  polynomials  Hili2_iN  (N  —  0, 1,  2, ... )  as  follows 


/  —  /(0)  ^clH  +  +  —ciijHij  +  •  •  •  +  ■j^O‘i1i2...iNHi1i2...iN  +  . . .  ^  ,  (112) 

where  ani2...iN  (N  =  0, 1,  2, ... )  are  tensorial  coefficients  that  depend  on  x  and  t.  Fur¬ 
thermore,  the  Maxwellian  distribution  function  /(0)  is  written  as 


<40, 


where  cn(^)  is  the  weight  function: 


<40 


1  r-€z/2 
(2tt)§ 


with  Q 


(113) 


(114) 


On  the  basis  of  the  weight  function,  the  tensorial  Hermite  polynomials  are  written  as 


-Hiii2...i;v(0 


(-1)*  dNuj 

u  dQ/JQz . . .  8Qn  ’ 


(115) 


being  orthogonal  with  respect  to  cu(£),  i.e., 


J  w(0  Hi  ■■Pm  mi  =  $MN  Aiijii2j2—iN3N-  (116) 

In  the  above  equation,  A iljli2j2...iNjN  represents  the  sum  A iljli2j2...iNjN  =  + 

(all  permutations  of  the  j's  indices)  and  5mn  is  Kronecker’s  symbol. 
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From  (115)  we  can  obtain  the  first  four  tensorial  Hermite  polynomials,  namely, 


m) 


l,  //,•(£)  -  o-  //p(C)  -  GA.;  -  *ij. 
Hijk(€)  £,i^j^k  (^i$jk  T  ^j$ik  T  ^k$ij)- 


(117) 

(118) 


In  terms  of  the  tensorial  Hermite  polynomials  we  can  write  the  thirteen  moments  of 
the  distribution  function  as 

e  =  /m(^“)2  fHd^  °  =  /m(^)  fHidt  (119) 

T-wejm{^y  <i2o) 

P(ij )  =  J  m  f  (h*]  ~  \Hrr$ij^  d£,  ft  =  I  y  f  (Hin  +  5Ri)  (121) 


We  recall  that  the  temperature  is  given  by  T  =  2mprr/3kg. 

The  coefficients  aj1j2...jiV  that  appear  in  the  distribution  function  (112)  are  determined 
from  the  definition  of  the  moments  of  the  distribution  function.  Furthermore,  for  a  thirteen 
moment  theory  the  distribution  function  (112)  is  written  as 

/  =  /<»>  (aH  +  aiH,  +  .  (122) 


where  a,  a j,  and  arri  represent  thirteen  scalar  coefficients  to  be  determined.  We  note 

that  the  term  a^kH^^  was  decomposed  according  to 


1 

6 


a(ijk )  +  g  (  drri^jk  T  ^rrj^ik  T  Q'rrk&ij) 


Hijk  Q^(ijk)  Hijk  T  ^  q  drrl Hss{ 5 


(123) 


and  the  part  associated  with  the  third  order  tensor  | a^j^Hijk  was  not  taken  into  account. 

If  we  insert  now  the  distribution  function  (122)  into  the  definitions  of  the  moments 
(119)  -  (121),  integrate  the  resulting  equations  and  use  the  orthogonality  properties  of 
the  tensorial  Hermite  polynomials  (116),  we  obtain  that 


a  =  1,  a%  —  0,  arr  =  0, 


„  _  PM  _  2g*  m_ 

a(d>  p  i  arrt  p  \j  kT' 


Hence,  the  distribution  function  for  the  thirteen  moments  (with  =  m/2kT)  becomes 


/  =  f,0>  { i  +  — 

e 


P(ij)C,Cj  +  ig.a  I  I3C!2  -  5 


(125) 


which  is  the  so-called  Grad’s  distribution  function. 

In  order  to  determine  the  non-convective  fluxes  pijk  and  q, ij  we  insert  Grad’s  distribu¬ 
tion  function  (125)  into  the  definitions  (110)i  and  (lll)i  and  obtain  by  integrating  the 
resulting  equations: 


Pijk  T  qjSik  T  qkSij)  * 


_  V  x  ,  7P 

%  2 e  Sii  2gm' 


(126) 
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For  the  calculation  of  the  production  terms  Pij  and  Qi  we  change  the  integration 
variables  (c,  Ci)  i — >  (g,  G)  so  that  in  these  new  variables  (110)2  and  (111)2  become 


p,  =  /  - 


2  L2 


,,  ~  9i9j )  -  Gi  (g'j  -  (jj)  -  Gj  {g[  -  </,-) 


ffidT, 


(127) 


m 


Qi~  /  2 


(9r9i  ~  9r9i )  ~  x  (  G2  +  -g 2  )  (r/(  —  gi)  —  GiGr  (g'r  —  gr) 


ffidT,  (128) 


where  the  product  of  Grad’s  distribution  functions  in  a  linearized  theory  reads 


//,  =  +  4Tpw 


+  -gkgi 


+ 


5  Q 


9j 


2Gj  (  /3G2  + 


o  +  tdgjGkgk 


The  integration  of  (127)  and  (128),  yields 
P{ij ) 


P  = 
1  v 


Tr 


Qi  =  ~, 

6Tr 


where 


m 


Tr  = 


5 


16n  V  vr kT  0(2’2)  ’ 


(129) 

(130) 

(131) 


represents  the  relaxation  time  of  the  pressure  deviator. 

Once  we  know  the  constitutive  equations  for  the  non-convective  fluxes  (126)  and  pro¬ 
duction  terms  (131)  as  functions  of  the  thirteen  basic  fields  g,  Vi,  and  q,,  a  system  of 
linearized  thirteen  scalar  field  equations  can  be  obtained  from  (106)  -  (109),  namely 


Vo 

+ 

d  Vi 

=  °, 

OXi 

(132) 

gVvi 

+ 

dp  dp{ij)  _ 
dxi  dxj 

(133) 

nkVT 

+ 

dqi  ,  diq  dvi 

dTi+pdTpp™dTi-0’ 

(134) 

Vp(ij) 

+ 

A9(i(i  ,  y9v(i  _  P(ij) 

5  dx^  y  dx^  Tr  ’ 

(135) 

Vqi 

+ 

p  dpgk)  5  k  dT  _  2  qi 
g  dxk  2  mP  dxi  3 rr 

(136) 

If  we  restrict  ourselves  to  a  five  field  theory  described  by  the  fields  of  mass  density, 
momentum  density  and  temperature,  the  corresponding  balance  equations  are  given  by 
(132),  (133)  and  (134).  In  this  case  the  pressure  deviator  and  the  heat  flux  vector  are  no 
longer  variables,  just  constitutive  quantities.  We  may  use  the  remaining  eight  equations 

(135)  and  (136)  in  order  to  obtain  the  constitutive  equations  for  the  pressure  deviator 
and  for  the  heat  flux  vector,  by  considering  only  the  equilibrium  values  of  the  pressure 
deviator  =  0  and  of  the  heat  flux  vector  g*  =  0  on  the  left-hand  sides  of  (135)  and 

(136) .  Hence,  only  the  underlined  terms  remain,  and  we  obtain 


P(ij) 


—2  g 


dv(j 
dxj)  ’ 


9i  =  -A 


dT 

dxi 


where  // 
where  A 


5  ImkT  1 

PTr  =  16  V  vr  IT2-2)’ 

15 k  _  75 k  ImkT  1 

4 mP  r  64 mV  vr  IT2’2) 


(137) 

(138) 


9-22 


RTO-EN-AVT-194 


6  CHAPMAN-ENSKOG-GRAD  COMBINED  METHOD 


Equations  (137)  and  (138)  represent  the  laws  of  Navier-Stokes  and  Fourier,  respectively, 
and  the  coefficients  of  shear  viscosity  /i  and  thermal  conductivity  A  correspond  to  the  first 
approximation  to  these  coefficients  found  by  applying  the  Chapman-Enskog  method  (see 
(68)). 

6  Chapman-Enskog-Grad  Combined  Method 

We  describe  now  another  method  to  obtain  the  constitutive  equations  for  the  pressure 
deviator  and  for  the  heat  flux  vector  which  combines  the  features  of  Chapman-Enskog 
and  Grad’s  methods.  In  this  method  neither  a  solution  of  the  integral  equation  is  needed 
-  as  in  the  Chapman-Enskog  method  -  nor  the  field  equations  for  the  moments  of  the 
distribution  function  are  used  -  as  in  Grad’s  method. 

To  begin  with  we  observe  that  in  the  representation  /  =  /(0^(1  +  0)  the  deviation  of  the 
Maxwellian  distribution  function  <f>  in  the  Chapman-Enskog  method  -  which  is  a  unknown 
quantity  -  can  be  written  in  terms  of  known  quantities.  These  known  quantities  can  be 
chosen  as  the  thirteen  scalar  fields  of  mass  density,  hydrodynamic  velocity,  temperature, 
pressure  deviator  and  heat  flux  vector  which  appear  in  Grad’s  distribution  function  (125), 
namely, 

9 S2  T  4  /  5\' 

</>  = P{lj)CiCj  +  -qiCi((3C2-- j  .  (139) 

In  the  framework  of  the  Chapman-Enskog-Grad  combined  method  we  insert  the  rep¬ 
resentation  (139)  into  the  non- homogeneous  integral  equation  of  the  Chapman-Enskog 
method  (34),  so  that  it  becomes  an  equation  for  the  determination  of  the  pressure  devia¬ 
tor  p^p  and  of  the  heat  flux  vector  qt.  This  equation  reads 

/<»>  { 1  (>  -§)<£  +  =  frm  mn] 

+  |^,I  (V'2  -  Q  .  (140) 

For  the  determination  of  the  pressure  deviator,  we  multiply  (140)  by  CqCp  and  inte¬ 
grate  the  resulting  equation  over  all  values  of  c,  yielding 

v  dvr  2/9^  P 

ePm  .r{,Ci>IlCkC,]dc-  (141) 

To  solve  the  above  equation  for  the  pressure  deviator  ,  we  introduce  the  integral 

Iim  =  I  CiGZ[GC,]dc  =  +  SuSik),  (142) 

which  was  written  in  terms  of  Kronecker’s  symbol  with  coefficients  A\  and  A2.  If  we 
contract  the  integral  Iijki  in  two  different  manners,  namely,  Irrss  =  9A±  +  6A2  and  Irsrs  = 
3Ai  +  12A2  we  get  that 

—  T 

oirsrs  J-rrss  r  r  (1  AQ\ 

t  {ij)kl  =  - ^ 
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We  use  now  the  relationship  (143)  to  write  (141)  as 

2-^  =  If-pw  [ C{rC,}X[C{rCt)]dc  = 


m  dx 


j) 


5  Q 


5  Q 


(144) 


where  /3 (°’°)  is  given  by  (65).  Hence,  it  follows  from  (144)  Navier-Stokes’  law 

dv(i  ,  5  ImkT  1 


P(ij)  2/r 


dxj)  ’ 


where 


li  = 


16 


7T 


(145) 


is  the  coefficient  of  shear  viscosity. 

We  determine  the  heat  flux  vector  by  following  the  same  methodology.  We  multiply 
(140)  by  C2Ci  and  the  integrate  of  the  resulting  equation  over  all  values  of  c,  resulting 


kp  d T 


qk  /  C2QX 


m2  dxi  5  q 

thanks  to  the  following  relationship 

iik  =  f  c2ca 


(3C2 


\)Ck 


dc  = 


15  g 


a(1,1)  ft, 


-  -  )  Ck 


dc  =  —d 


ik* 


Hence,  Fourier’s  law  follows  from  (146) 

where 


x  dT 
Qi  =  -Aw—, 
OXi 


75  k  mkT  1 
64  m  V  vr  fA2’2)  ’ 


(146) 


(147) 


(148) 


is  the  coefficient  of  thermal  conductivity.  Above,  we  have  used  the  relationship  cA1’1)  = 
^(°>°)  which  was  obtained  in  the  Chapman- Enskog  method. 

We  can  extended  this  method  to  obtain  the  successive  approximations  to  the  transport 
coefficients,  by  introducing  traceless  second  order  tensors  pfl  and  vectors  qf'1  defined  by 


(s)  (2s)!! 

4>  = 15- 


(2s +  5)!! 


mS^C^fdc,  #  = 


15(2s)!! 


4/3  (2s  +  3)!! 


mSPCifdc,  (149) 
2 


where  =  pqp  represents  the  pressure  deviator  and  qf]  =  q,  the  heat  flux  vector, 

while  the  the  higher-order  moments  pf^ (s  >  0)  and  qf'1  (s  >  1)  do  not  have  standard 
designations. 

For  the  moments  of  the  distribution  function  characterized  by  g,  T,  pff  and 
the  distribution  function  reads 


/  =  fm  \ 1  +  v  E  -  f)  E  sf  cd* 

^  s=0  2  ^  s=0  2 


(150) 


so  that  (140)  becomes 


1 


/<"’  1  -4sfc^ + \  -  4  E  AD 


dT 


T 


dxk 


do) 

’5 

2 


dv 


(k 


dxr 


2f2 


» 


s= 0 


o(s)  s~i  n 
>->5  U-fcO-i 


5g 


E«A 


s=0 


q(s)s~i 
*->3  C'k 
2 


(151) 
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The  multiplication  of  (151)  by  /3S^CuCj\  and  integration  of  the  resulting  equation 

2 

leads  to 


5  O'  5(0,r)dv(i 
2  m/3 3  dxj) 


s=0 


(152) 


If  we  solve  the  above  equation  for  the  pressure  deviator  ppp  it  follows  the  Navier-Stokes 
law 

PM  =  -2  H  where  H  =  (/5-1)(0’°)  (153) 

is  the  coefficient  of  shear  viscosity  in  the  pth  successive  approximation.  In  (153)2  (/3“1)<'°'°') 
represents  the  first  row  and  first  column  of  the  inverse  of  the  p  x  p  matrix  j3^r,s\ 

If  we  multiply  (151)  by  and  integrate  the  resulting  equation  we  get 


75  n2k  (l  r)  d T 
16  (3  dxi 


oo 


s=0 


which  can  be  solved  for  the  heat  flux  vector  (p,  yielding  Fourier’s  law: 


Qi  =  ~  M 


dT 

p  d  Xi 


where 


75  n2k 

16~/T 


-n(i-i) 


(«") 


(154) 


(155) 


The  coefficient  of  thermal  conductivity  [A]p  in  the  pth  successive  approximation  is  given 
in  terms  of  the  element  (a^1)^1,1^  of  the  inverse  matrix. 

The  successive  approximations  of  the  transport  coefficients  of  shear  viscosity  (153)2 
and  thermal  conductivity  (155)2  are  the  same  as  those  which  follow  from  the  Chapman- 
Enskog  method. 


7  Fourteen  Moment  Theory  for  Granular  Gases 

We  may  characterize  a  macroscopic  state  of  the  granular  gas  by  the  fourteen  fields  of  mass 
density  g,  hydrodynamic  velocity  vu  pressure  tensor  p,jj ,  heat  flux  vector  q \  and  contracted 
fourth  moment  pajj  which  are  defined  by 


6=  mfdc,  QVi=  mcifdc,  ptj  =  mQCjf  dc, 


m 


Qi  =  /  y C  Cif  dc>  Piijj  =  /  mC  f  (lc- 


(156) 

(157) 


The  balance  equations  for  the  fourteen  fields  are  obtained  from  the  Boltzmann  equation 
(78)  by  multiplying  it  successively  by  m,  mci,  mCiCj,  mC2Cij 2  and  mCA  and  integrating 
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the  resulting  equations.  Hence,  it  follows 

dp  dpvi  = 
dt  dxi  ~  ’ 

dgi’i  d  ,  . 

-5r  +  s-(^+p«)  =  o, 


(158) 

(159) 


7T  +  a7  (p«‘ + p«"‘> +  P“I|  +  p«t)  =  p'”  (160) 

dqi  d  dvj  dvi  Pkidpkj  lprrdpij 

7  +  7“  +p«‘§7  +  ®a^  -  7  "  27  &;  =  0i’  (161) 


<9  .  dvl  8  (9p*fc 

A  TUT  yPiiijk  A  PiijjVj)  A  ^1ik~aZ  “'h  o  =  A*- 


dt  dx 


dxk  q  dxk 


(162) 


In  the  above  equations  the  new  moments  of  the  distribution  function  are  defined  by 

Pijk  =  I  mCiCjCkfdc,  qi:j  =  J  ™C2CiCjfdc,  piijjk  =  J  mC4Ckfdc,  (163) 
where  qa  =  Pujj/ 2.  Furthermore,  the  production  terms  read 

Pij  =  \  J m  (cl' a j'  +  C'C'  -  C\C)  -  C'C'  )  /i/d2  (g  •  k)  dk dc,  dc,  (164) 

Qi  =  \  [  j  (c?Cty  +  G'2G'  -  Cl  Cl  -  G2G,:)  fifd2  (g  •  k)  dk  dCl  dc,  (165) 


P=-  m  (C4  +  G'4  -  G4  -  G4)  /i/d2  (g  •  k)  dk  dc,  dc. 


(166) 


We  can  decompose  the  balance  equation  for  the  pressure  tensor  (160)  in  its  trace  and 
traceless  parts  by  introducing  the  traceless  tensors  pqj)  and  pq-jk)  defined  by 


Pij  P(ij)  A  pSij,  Pijk  P(ijk)  A  ( Qidjk  A  qj&ik  A  qk$ij)  j  (167) 

where  p  =  nkT  is  the  hydrostatic  pressure.  Hence,  we  obtain  from  (160)  the  equations: 


dT  ,  dT  2  ( dqi  dvi  diq  \  _ 

s^  +  ”‘ay  +  3 ^(ay+pay+PW>syj+TC-0’ 

<9p(P)  d  (  x  ,  Ac,;  2  chy 

WW  +  7G7  W)  +P(u>wO 


dt  dxk 
Adgg  dv{i  _ 

5dxj )  J  dx^  ^ ' 


(168) 


(169) 


Equation  (168)  is  the  balance  equation  for  the  temperature  with  /  denoting  the  cooling 
rate  (81),  while  (169)  is  the  balance  equation  for  the  pressure  deviator. 

We  represent  the  non-eqnilibrium  distribution  function  for  the  fourteen  moments  as 


/  —  /m  (a  A  a,iCi  +  ciijCiCj  +  biC2Ci  +  bC4)  , 


(170) 


where  the  fourteen  coefficients  a,  a*,  ai?-,  and  b  do  not  depend  on  the  peculiar  velocity. 
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At  equilibrium  the  contracted  fourth  moment  reads 

Pan  =  J  mClfM  dc=lbe(^)  , 


(171) 


and  it  is  convenient  to  introduce  a  dimensionless  non-equilibrium  part  of  pujj,  denoted  by 
A  and  defined  by 


A  —  — J—  f— V 
15p  \kT J 


(o) ' 

Piijj  Piijj 


-  i mcv  - Mdc-  (m) 


The  fourteen  coefficients  a,  a*,  a^-,  b  are  determined  from  the  definition  of  the  basic 
holds  (156),  (157)  and  (172)  together  with  the  non-equilibrium  distribution  function  (170), 
yielding 


2/9 


8/32 


/  =  M  1  +  —p^CrC,  +  —  q^  (  f3C2 


15  5  (3C2  (32C4\ 


+ 


J 


A 


(173) 


We  can  now  determine  the  constitutive  equations  for  the  moments  of  the  distribution 
function  by  inserting  (173)  into  the  expressions  (163)  and  subsequent  integration  of  the 
resulting  equations.  Hence,  it  follows 

kT  5  f kT\ 2  7  kT 

P{ijk)  =  0,  pnjjk  =  28—qk,  Qij  ~  2  ^  \~rn  J  ^  ^  ^  2 

The  constitutive  equations  for  the  cooling  rate  (  and  for  the  production  terms  P(ij),  Qi ,  P 
are  obtained  in  the  same  manner  and  it  follows  from  (81),  (164),  (165)  and  (166)  through 
integration  that 


c 


1 

3  r 


1 

5  r 


(1  +  ck)  (3 


a)P{ij)i 


Qi  =  +  a)  t49  _  33“] 

(1-o2)(2 a2  +  9)v/L{i  + 


30a2(l  -  a)  +  271  -  207a  A  ) 
(2ck2  +  9)(1  —  a)  16  J  ' 


(175) 

(176) 

(177) 


We  have  introduced  in  the  above  equations  a  mean  free  time  in  terms  of  a  reference 
temperature  To,  namely,  r  =  ^/^r/4nd2.  Furthermore,  we  have  considered  in  the  above 

expressions  only  linear  terms  in  A.  For  the  elastic  case  a  =  1  and  the  production  terms 
(175)  -  (177)  reduce  to 


C  =  0,  P{ij) 


4  Ft 

~5r  V  T^P{iJh 


Qi 


8  f Y 
15r  V  To^ 


P 


8p  ( kT\ 
t  \m  J 


T 

T0 


A.  (178) 


From  the  knowledge  of  the  constitutive  equations  in  terms  of  the  basic  fields,  we  get 
-  by  inserting  (174)  -  (177)  into  the  balance  equations  (158),  (159),  (161),  (162),  (168) 
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and  (169)  -  the  following  linearized  system  of  fourteen  field  equations  for  g,  u,,  T,  P(ij)„ 
Qi  and  A: 


_  d  Vi 

Vq+  Qt—  =  0, 

OXi 


eVv,  +  %L  +  9pM  =  o, 

OXi  OXi 


VT 


dq.i  dvi\  T  j~T  2. 

smtV&W&pWvU  1 


1  +  A 
16 


=  0, 


4  dqu  d v/i  1  [F .  .  . 

Vpm  +  ld^~)+  2po. T  =  “57  V  A1  +  “)( 


a 


)P(ij)i 


m  9x,- 


2m 


<9A  dn  \  5 k2Tn  d T 

n 7w-  +  A—  + 


<9au 


dxi) 


2m  dxi 


60 r  V  To 
A;T\  2 
m  / 


—  (1  +  a)  [49  -  33a]  qu 
8  kT  dq 


m  dxi 


t  \m  ) 


(179) 

(180) 
(181) 


(182) 


15e|—  ]  UA+—  (l  +  a)A/l^(l-2a2)(l-a) 

-*0 


A 


[81  -  17ck  +  30a2(l  —  a)]  — 

16  1 


(183) 


In  (183)  we  have  used  (179)  and  (180)  in  order  to  eliminate  the  time  derivatives  of  o  and 
T  and  neglected  products  of  A. 

Now  we  search  for  spatially  homogeneous  solutions  of  the  fourteen  held  equations 
where  the  holds  depend  only  on  time.  The  held  equations  (179)  -  (183)  in  this  case  read 

(184) 

(185) 

(186) 
(187) 

[81  -  17a  +  30a2(l  -  a)Al .  (188) 

16  J 

Above,  the  following  dimensionless  quantities  were  introduced:  time  A  =  t/r,  temperature 
T*  =  T/T0,  pressure  deviator  p*^  and  heat  hux  vector  q*. 

We  may  infer  from  (184)  that  the  mass  density  and  the  velocity  helds  remain  constant 
in  time,  while  (185)  -  (188)  compose  a  system  of  coupled  differential  equations  for  the 
determination  of  the  temperature  T*,  pressure  deviator  pF,  heat  hux  vector  q*  and 
fourth  moment  A.  Since  this  system  of  equations  is  non-linear,  it  was  solved  numerically 
by  considering  the  initial  conditions  T*(0)  =  1,  pW( 0)  =  1,  q*( 0)  =  1  and  A(0)  =  1.  The 
left  frame  of  Figure  3  shows  the  time  decay  of  the  temperature  (solid  line)  in  comparison 
with  Haff’s  law  (dashed  line)  which  is  the  solution  of  (185)  when  A  =  constant  (see  (190) 
below).  We  may  conclude  that  the  temperature  decay  T*  follows  closely  Haff’s  law  and  by 


de_  =  0  dVi  _  Q 

dA  dA 


3 

(IT  T2 

U±*  I  ±*  7i  2 

77  +  T(1-“ 


i  +  — 
16 


=  0, 


dp 


{ij) 


(1  +  a)(3  —  a) 


Vr*p*{ij) , 


dA 

dq*  _  (1  +  a) 

dA  60 


Vt,  [49  -  33a]  q*, 

f 
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increasing  the  restitution  coefficient  the  time  decay  of  the  temperature  is  less  accentuated. 
Furthermore,  from  the  right  frame  of  the  Figure  3  we  may  observe  that  the  pressure 
deviator,  the  heat  flux  vector  and  the  dimensionless  fourth  moment  decay  also  with  time 
and  the  trend  to  equilibrium  is  more  pronounced  for  the  pressure  deviator  followed  by 
the  heat  flux  vector  and  the  dimensionless  fourth  moment.  The  curves  in  the  right  frame 
were  obtained  for  the  restitution  coefficient  a  =  0.75.  By  increasing  the  value  of  the 
restitution  coefficient  the  time  decay  of  the  pressure  deviator  and  dimensionless  fourth 
moment  are  more  accentuated  but  the  one  of  the  heat  flux  vector  is  less  pronounced,  since 
it  is  connected  with  the  transport  of  energy  which  is  directly  affected  by  the  inelasticity. 


Figure  3:  Left  frame:  time  decay  of  the  temperature.  Right  frame:  time  decay  of  the 
pressure  deviator,  heat  flux  vector  and  fourth  moment. 


From  the  analysis  of  the  differential  equations  (186)  we  note  that  the  pressure  deviator 
and  the  heat  flux  vector  do  not  evolve  with  respect  to  time  when  the  initial  conditions  for 
these  fields  vanish.  However,  a  vanishing  initial  condition  for  A  implies  from  (188)  that  it 
could  evolve  with  time.  An  interesting  case  is  the  one  where  the  fourth  moment  remains 
constant  in  time  and  which  is  expressed  by  the  condition  that  the  right-hand  side  of  (188) 
must  vanish,  i.e., 


16(1  —  a)(l  —  2a2) 
30a2(l  —  a)  +  81  —  17a 


(189) 


The  above  expression  for  A  is  the  same  as  that  found  for  the  coefficient  a-2  which  follows 
from  the  Chapman- Enskog  method  (see  (96)).  We  note  that  A  vanishes  in  the  elastic  case 
where  a  —  1. 

If  we  substitute  the  expression  of  A  given  by  (189)  into  (186)  and  the  integrate  the 
resulting  equation,  we  get  Half’s  law 


T*(t) 


1  + 


6 


1  + 


3(l-a)(l-2a2) 

81-17a+30a2(l-a) 


(190) 


By  considering  the  elastic  case  a  —  1,  the  temperature  T*  remains  constant  in  time 
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and 


(191) 


i.e,  both  decay  exponentially  with  time,  which  is  a  well  known  result. 

In  the  case  of  a  five  field  theory  described  by  the  fields  of  mass  density,  velocity 
and  temperature,  the  pressure  deviator  and  the  heat  flux  vector  are  no  longer  fields  but 
constitutive  quantities.  The  constitutive  equations  for  these  quantities  are  obtained  from 
the  eight  remaining  field  equations  (181)  and  (182)  by  inserting  the  equilibrium  values 
P(ij)  =  0  and  Qt  —  0  on  the  left-hand  side  of  these  equations,  yielding  the  laws  of  Navier- 
Stokes  and  Fourier 


P(ij )  2/X 


dv{i 


dx 


j) 


.  d  T  d  n 
Qi  =  “Ax - v— , 

OX;  OX; 


(192) 


respectively.  The  transport  coefficients  of  shear  viscosity  n,  thermal  conductivity  A  and 
the  coefficient  d  are  given  by 


5  ImkT  1  75  k  ImkT  1 

^  4d2  V  7 r  (1  +  ct)  (3  —  a)  ’  2d2mV  n  (1  +  a)  (49  —  33a:)  ’ 

n  _  75 kT  ImkT _ 16(1  -  a)(l  -  2a2) _ 

2nd2mV  vr  (1  +  a)  (49  -  33a)  [30a2  (1  -  a)  +  81  -  17a] '  1  ’ 

We  note  that  if  we  consider  a  thirteen  moment  theory  the  heat  flux  vector  is  proportional 
only  to  the  temperature  gradient  and  the  dependence  on  the  particle  number  gradient 
does  not  show  up.  This  fact  can  be  understood  by  observing  that  this  dependence  comes 
out  from  the  underlined  term  in  (182)  which  depends  exclusively  on  the  dimensionless 
fourth  moment.  Furthermore,  the  transport  coefficients  (193)  and  (194)  are  not  the  same 
as  those  obtained  from  the  Chapman- Enskog  method  (104)  and  (105),  and  it  seems  that 
the  equivalence  of  the  two  methods  can  be  attained  if  we  consider  more  moments  of  the 
distribution  function  in  Grad’s  moment  method. 
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